#include <CGraph.h>
#include <fstream>
#include <string>
#include <iostream>
#include <CVector.h>

void CGraph::FFT() {
  double *in;
    fftw_complex* out;
    fftw_plan p;

    in = (double*) fftw_malloc(sizeof(double) * Bins);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Bins);
    
   
    for (int i=0;i<Bins;i++) {
      in[i] = Val[i];
      out[i][0] = 0.0;
      out[i][1] = 0.0;
    }

    p = fftw_plan_dft_r2c_1d(Bins, in, out,FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
   
    for (int i=0;i<Bins;i++) 
      ValFFT[i] = complex<double>(out[i][0],out[i][1]);

   
    fftw_free(in); 
    fftw_free(out);
  }

void CGraph::FFTInv() {
    double *out;
    fftw_complex* in;
    fftw_plan p;
    
    out = (double*) fftw_malloc(sizeof(double) * Bins);
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Bins);
    
    
    for (int i=0;i<Bins;i++) {
      in[i][0] = ValFFT[i].real();
      in[i][1] = ValFFT[i].imag();
      out[i] = 0.0;
    }

    p = fftw_plan_dft_c2r_1d(Bins, in, out,FFTW_ESTIMATE);

    fftw_execute(p); /* repeat as needed */
    fftw_destroy_plan(p);

    for (int i=0;i<Bins;i++) {
      Val[i] = out[i] / (double)Bins;
    }
    
   
    fftw_free(in); 
    fftw_free(out);
  }


void CGraph::SaveBin(string Filename) {
    ofstream f(Filename.c_str(), ios::out | ios::binary);
    f.write((char *)&Bins, sizeof(unsigned int));
    //    f.write((char *)Val, sizeof(double)*Bins);
    f.close();

  }

void CGraph::LoadBin(string Filename) {
    ifstream f(Filename.c_str(), ios::in | ios::binary);
    f.read((char *)&Bins, sizeof(unsigned int));
    Initialize(Bins);
    //f.read((char *)Val, sizeof(double)*Bins);
    f.close();
  }


void CGraph::cutModes(int from, int to) {
  for (int i=from; i<to; i++) {
    ValFFT[i] = complex<double>(0,0);
  }
}

void CGraph::Mean() {
  mean = 0;
  for (int i=0;i<Bins;i++) {
    mean+=Val[i];
  }
  mean/=(double)Bins;

}

void CGraph::Std() {
  std = 0;
  Mean();

  for (int i=0;i<Bins;i++) {
    double v= Val[i] - mean;
    std += v*v;
  }
  std/=(double)Bins;
  std = sqrt(std);

}

void CGraph::Scale(double v) {
  for (int i=0;i<Bins;i++) {
    Val[i] *= v;
  }

}


void CGraph::scaleX(double v) {
  for (int i=0;i<Bins;i++) {
    IndexScaled[i] = Index[i]*v;
  }

}


void CGraph::Gaussian(double B) {
  int x[2] = {0,1};

  double sum = 0;
  double scale = (Index[Bins-1] - Index[0])/(double)Bins;
  //cout << scale << endl;
  for (int i=0;i<Bins; i++) {
    double dist = 1E10;
    double val = 0;
    for (int l=0;l<1;l++) {
      double d = scale*(0 - i);
      //d/=Bins;
      dist = min(dist,d);
      val += exp(-dist*dist/(2.0*B*B));
      
    }
    
    sum+=val;
    Val[i] = val;
  }
  for (int i=0;i<Bins;i++)
    Val[i]/=sum;

}


void CGraph::Smooth(double B) {
  
  CGraph gauss;
  gauss.Copy(*this);
  gauss.Gaussian(B);
  gauss.SaveText("gauss.txt");
  gauss.FFT();
  FFT();
  MulFFT(gauss);
  FFTInv();

  

}


void CGraph::dampNoise(double A, double B) {
  for (int i=0;i<Bins;i++) {

    double scale = exp(-A*(double)i);
    Val[i] *= scale;
  }
  
}

void CGraph::LoadText(string Filename) {
  ifstream f(Filename.c_str(), ios::in);
  if (!f.good())
    throw string("File " + Filename + " not found");
  char inp[2000];
  Bins = 0;
   while(!f.eof()) {
      f.getline(inp,2000);
      Bins++;
     }
   Bins--;
  f.close();
  Initialize(Bins);
  f.open(Filename.c_str(), ios::in);
  int i= 0;
   while(!f.eof()) {
      f.getline(inp,2000);
      vector<string> tok;
      CUtil::Tokenize(inp, tok," ");
      if (tok.size()==2) {
	Index[i] = atof(CUtil::trim(tok[0]).c_str());
	IndexScaled[i] = Index[i];
	Val[i] = atof(CUtil::trim(tok[1]).c_str());
	i++;
      }
     }

  f.close();
}



void CGraph::SaveText(string Filename) {
  ofstream f(Filename.c_str(), ios::out);
  for (int i=0;i<Bins;i++) {
    f << IndexScaled[i] << " " << Val[i] << endl;
  }
  f.close();
}

void CGraph::SaveTextFFT(string Filename) {
  ofstream f(Filename.c_str(), ios::out);
  for (int i=0;i<Bins;i++) {
    f << Index[i] << " " << norm(ValFFT[i]) << endl;
  }
  f.close();
}


void CGraph::RenderGauss(double position, double sigma, double amp) {
  for (int i=0;i<Bins;i++) {
    double p = position - Index[i];
    p = pow(p/sigma, 10);
    double val = amp*exp(-p);
    Val[i] = max(Val[i], val);
  }
}



float CGraph::getValAtScaledIndex(double& index) {
  for (int i=0;i<Bins-1;i++) {
    if (index>=IndexScaled[i] && index<IndexScaled[i+1])
      return Val[i];
  }
  return NONE;
}

double CGraph::ChiSQ(CGraph& temp, CGraph& two) {
  double chisq = 0;
  for (int i=0;i<temp.Bins;i++) {
    double val = two.getValAtScaledIndex(temp.Index[i]);
    //double val = two.Val[i];
    if (val!=NONE) {
      chisq += abs(pow((temp.Val[i] - val), 2.0));
    }
  }
  return chisq;
}



int CGraph::getIndexAtScaledIndex(double& index) {
  for (int i=0;i<Bins-1;i++) {
    if (index>=IndexScaled[i] && index<IndexScaled[i+1])
      return i;
  }
  return -1;
}



void CGraph::MulWithIndex(CGraph& o) {
    for (int i=0;i<Bins;i++) {
      int j = o.getIndexAtScaledIndex(IndexScaled[i]);
      if (j!=-1) {
	Val[i] *= o.Val[j];
      }
    }
  }


void CGraph::calculateStatistics() {
  double top = -1;
  int idx;

  for (int i=0;i<Bins;i++) {
    if (Val[i]>top) {
      top = Val[i];
      idx = i;
    }
  }
  mean = Index[idx];

  double val = Val[idx];
  while (Val[idx]>0.5)
    idx--;

  double pos = Index[idx];
  std = abs(mean - pos);

}

